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We study by mean-field analysis and stochastic simulations chemical models for genetic toggle 
switches formed from pairs of genes that mutually repress each other. In order to determine the 
stability of the genetic switches, we make a connection with reactive flux theory and transition state 
theory. The switch stability is characterised by a well defined lifetime t. We find that r grows 
exponentially with the mean number A'^ of transcription factor molecules involved in the switching. 
In the regime accessible to direct numerical simulations, the growth law is well characterised by 
r ~ N°'exp{bN), where a and 6 are parameters. The switch stability is decreased by phenomena 
that increase the noise in gene expression, such as the production of multiple copies of a protein from 
a single mRNA transcript (shot noise), and fiuctuations in the number of proteins produced per 
transcript. However, robustness against biochemical noise can be drastically enhanced by arranging 
the transcription factor binding domains on the DNA such that competing transcription factors 
mutually exclude each other on the DNA. We also elucidate the origin of the enhanced stability 
of the exclusive switch with respect to that of the general switch: while the kinetic prefactor is 
roughly the same for both switches, the 'barrier' for flipping the switch is significantly higher for 
the exclusive switch than for the general switch. 

PACS numbers: 05.40.-a; 87.16.Yc 



I. INTRODUCTION 

In an organism, genes can be turned on or off by the 
binding of proteins to regulatory sites on the DNA in the 
vicinity of the starting point for transcription I]. The 
proteins are known as transcription factors and the DNA 
binding sites are known as operators. The process is an 
example of gene regulation. Transcription factors can 
turn genes off by stereochemical blockage of the binding 
of RNA polymerase (RNAp), or they can turn genes on 
by co-operative binding (recruitment) of RNAp 

Since transcription factors are proteins, they are coded 
for elsewhere on the genome. This means that transcrip- 
tion factors can regulate the production of other tran- 
scription factors, or indeed can regulate their own pro- 
duction. A highly complex network of biochemical re- 
actions can be built up, capable, in principle, of solv- 
ing arbitrarily complex computational problems 0, Q- 
The network is interfaced to the outside world by, e.g., 
signalling cascades from receptor proteins or by spe- 
cific interactions between transcription factors and small 
molecules such as metabolites. In this way, even rela- 
tively simple organisms such as E. coli can perform fairly 
complex computations such as integrating different envi- 
ronmental signals. In higher organisms, gene regulation 
networks lie at the heart of cell differentiation and devel- 
opmental pathways. 

Genetic toggle switches are an informative paradigm in 
this context j^, , 5j . They are regulatory constructs which 
select between two possible stable states, representing for 
example differentiation between two developmental path- 



ways. Perhaps the simplest kind of switch is one that is 
constructed from a pair of genes that mutually repress 
each other, as indicated in Fig.^ The switch of A-phage 
in E. coli is a naturally occuring example, which has been 
studied in much detail 4, 6, 7, 8J. Another well-known 
example is that of the human herpesvirus 3 (chickenpox 
or varicella-zoster virus) , which has a pathogenesis that is 
remarkably similar to A-phage; the virus lies dormant af- 
ter the initial infection, but can be triggered to re-emerge 
much later as shingles (herpes zoster). Synthetic toggle 
switches have also been constructed in vivo 

Genetic switches are usually flipped by external sig- 
nals. In the A-phage , for example, the switch is initially 
in the dormant (lysogenic) state but can be flipped into 
the active (lytic) state by the presence of the bacterial 
protein RecA. Such an induction event occurs when the 
cell starts to produce RecA to repair DNA damage as a 
result of, e.g., a burst of UV light. 

Importantly, genetic switches are often so stable that 
they remain in one state until the external trigger flips 
the switch. In wild-type A-phage spontaneous flips are 
extremely rare, occurring at a rate of perhaps as low as 
one spontaneous flip in 10"'^^ s -'^^^ important question, 
therefore, is: what are the design principles that allow the 
switch to attain such extreme stability in the presence 
of fluctuations and biochemical noise? This question is 
particularly relevant, because detailed modelling has sug- 
gested that the stability of A-phage cannot be explained 
on the basis of our current understanding 'sj . 

We have recently shown that mutual exclusion of com- 
peting transcription factors can drastically enhance the 
stability of genetic switches H^. Transcription factors 
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FIG. 1: (a) A genetic toggle switch can be formed from a pair 
of genes that mutually repress each other's expression, (b) If 
the genes are transcribed on opposite strands of the DNA, the 
upstream regulatory domains can overlap and interfere with 
one another. As a result, competing regulatory molecules 
mutually exclude each other. 



can mutually exclude each other, if their operator sites 
(partially) overlap (see Fig. ^) . It appears that nature 
has exploited the functional benefit of this spatial ar- 
rangement of operators, since a recent statistical analy- 
sis has revealed that this network motif of 'overlapping 
operons' is over- represented in the bacterium E. coli jTll |. 

In this paper, we extend the analysis of Ref. ^3 in sev- 
eral ways In the next section, we first describe the 
set of chemical reactions by which we model the genetic 
toggle switches (see Fig. 0. We then turn to a mean- 
field (chemical rate equation) analysis of the appearance 
of switching behaviour. This mean-field analysis shows 
that the region of bistability is significantly larger for the 
exclusive switch than for the general switch. In order 
to determine the life times of the stable steady states of 
the switches, we have performed direct stochastic simu- 
lations; these are described in section Hvl We then con- 
sider a number of refinements of the basic switch model. 
In particular, we explicitly take into account transcrip- 
tion and translation. We show that both 'shot noise' 
and fluctuations in the number of proteins produced per 
mRNA transcript markedly decrease the switch stability. 
However, we again find that the exclusive switch is still 
orders of magnitude more stable than the general switch. 
In this paper, we also elucidate the origin of the enhanced 
stability of the exclusive switch with respect to that of 
the general switch: while the kinetic prefactor is roughly 
the same for both switches, the 'free-energy' barrier for 
flipping is significantly higher for the exclusive switch 
than for the general switch This underscores our 

earlier observation that mutual exclusion can drastically 
enhance the robustness of genetic switches to biochem- 
ical noise ^3- I* also strongly supports the conjecture 
that regulatory control is one of the main evolutionary 
driving forces that shape the organisation of genes and 
operons on the genome 11]. 



II. MODEL SPECIFICATION 

Our starting point is a set of chemical reactions which 
model the processes involved in the toggle switches shown 
in Fig. n As chemical species, we introduce a pair of 
transcription factors (TFs) which can exist as monomers, 
A and B, or multimers, A„ and B^. The multimers are 
responsible for regulating gene expression and are allowed 
to bind to the genome at an operator site. Multimers 
are introduced in order to get sufficient co-operativity in 
the binding isotherms to make a genetic switch |l,l| . as 
described in more detail below. The state of the operator 
is represented by O, 0A„, etc, depending on the state of 
binding of the multimers. The chemical reactions are 

nA ^ A„, 77iB ^ Bm, (fcf , /cb) (la) 

+ A„^OA„, + B„^OB„„ (fco„,fcoff) (lb) 

0A„ + B„ I 0B™ + A„ ^ OA„B„, (fco„,fcoff) (Ic) 

0|OA„-^A, 0|OB,„-^B, (fcA),(fcB) (Id) 

A|B^0, (a^a),(mb) (le) 

Here we have adopted a condensed notation in which '|' 
indicates alternative sets of reactants and '^' indicates 
that the reactants are not destroyed by the reaction. For 
example the first of Eqs. I|ld|l . O | 0A„ ^ A, represents 
two reactions O O + A and OA^ 0A„ -I- A. 

The reactions in Eqs. ^ account for, respectively, the 
formation of multimers, the binding of TF multimers to 
the operator (Eqs. Ijlb|l and lfTc|) '). the expression of TF 
monomers, and the degradation of TF monomers. Re- 
pression of gene expression is implicit in Eqs. (|ld|) . thus 
A is expressed if and only if B^ is not bound, etc. Re- 
action rates are as indicated, and we define equilibrium 
constants for multimerisation, Kd ~ fcf/fcb, and operator 
binding, K^^ = Kn/koS- 

Whilst detailed and biologically faithful models can 
be constructed as has been done for the A-phage switch 
1^ 13 , the above model is intentionally as simple as pos- 
sible. We believe that such an approach is as important 
as detailed biological modelling in elucidating the basic 
physical principles behind switches. Thus, for example, 
we have condensed the details of transcription and trans- 
lation into a single reaction step in Eqs. Ijld|l . governed by 
rate coefficients and fee • In a later section we explore 
the possibility of refining the model at this point. 

It should be noted, however, that the design of the 
network has to obey certain constraints. In particular, 
the TFs must bind co-operatively to the DNA in order 
to make a working switch In the present model, 

co-operativity is introduced through the binding of TF 
multimers rather than monomers. 

In our model the operator is in one of four states 
{0, 0A„, OBm, OA„Bm}- We now include the effect of 
interference between the upstream regulatory domains 
by disallowing some of these states. This is in the spirit 
of simplicity, strictly speaking the effect is to modify the 
probabilities of the states. Since the empty operator state 
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TABLE I: Distinct possibilities for the subsets of operator 
states for our switch modeh '0' indicates the state is disal- 
lowed, '1' indicates it is allowed. 



operator states 


O 


OA„ 




OA^iB^n 


case 




a 


P 


7 


general 


1 


1 


1 


1 


exclusive 


1 


1 


1 





partially co-operative 


1 


1 





1 


totally co-operative 


1 








1 



is always a possibility and both A„ and should be al- 
lowed to bind otherwise they would not be TFs, it turns 
out that there are only five possibilities, two of which are 
related by symmetry. The four distinct cases are shown 
in Table and are implemented by excluding some of 
the reactions in Eqs. (|lb|) and lfTc|) . For example, the 
exclusive switch is obtained by discarding the reactions 
in Eqs. lfTc|l thereby removing the state OA„Bm. 



III. MEAN-FIELD ANALYSIS 

We first analyse the behaviour of Eqs. Q using chem- 
ical rate equations. This plays the role of a mean field 
theory for this problem since chemical rate equations de- 
scribe the temporal evolution of the mean concentrations 
of molecules. Switching behaviour corresponds to the 
appearance of two distinct stable fixed points (attrac- 
tors) in the space of concentration variables. For the 
general switch, the problem has been analysed in detail 
by Cherry and Adler For the exclusive switch, a 

specific example has been studied by Kepler and Elston 
[iJI . Our approach is a generalisation of the analysis of 
Cherry and Adler. 

Let us consider what happens for A. Write for the 
number of A monomers in the cell etc, and let the cell 
volume he Vc- In a steady state, the multimcrisation 
reaction nA ^ A„ is in equilibrium and 



V V, 



(2) 



Similarly the binding reaction O -I- A„ ^ 0A„ (if it is 
allowed) is in equilibrium and 



( 



noA„ 



V Vc 



(3) 



Now, no is the probability that the operator is in state 
O times the number of copies of the genome in the cell 
(usually assumed to be one), etc. Therefore 



Prob(OA„) ^ K^IU 
Prob(O) 



(4) 



where we have introduced reduced concentration 

variable, equal to {l/Vc){KbKdY/"' times the number 



of monomers of A in the cell. Similarly we introduce 
a reduced concentration variable y for the number of 
monomers of B. 

^From the totality of binding equilibria, we sur- 
mise that the probabilities of the operator states 
{O, 0A„, 0B,„, OA„B„,.} are in the ratio 1 : as" : /3y™ : 
ry^riym .^^j^gj-g ^ffQ iiave covered off all the switch con- 
struction possibilities by introducing a set of coefficients 
(a, (3, 7) which take the values zero or one according to 
Table m The probability of the operator being in a state 
where A is expressed is therefore 



1 + ax" + /3y™ -I- 7a;"?/™ 
and of being in a state where B is expressed is 

1 -h/3y" 



9{x,y) 



1 + ax" + Py"^ + 7x"?/" 



(5) 



(6) 



For A, expression occurs at a rate fcA and degradation 
at a rate /^a x tia. It is convenient to define a reduced 
degradation rate 



MA 



Vc 



MA 



(i^bi^d)!/" kA ' 



(7) 



If there is more than one copy of the genome in the cell, 
we should additionally divide by the number of copies of 
the genome. A similar reduced degradation rate mb is 
defined for B. 

At a fixed point (a steady state) , the rate of expression 
is equal to the rate of degradation (for example kxf = 
Ma'T'a), for both TFs. Using the reduced concentrations 
and degradation rates defined above, the steady states 
are defined by 



f{x,y) = tt.AX, g(x,y) = MBy- 



(8) 



We see that the steady states only depend on /xa and mb ■ 
The simplest way to proceed now is to consider the 
fixed points for the dynamical system 



fix,y)-ilAX, ij ^ g{x,y) - fi^y. 



(9) 



This is only an approximation to the full dynamics, since 
it assumes that TF binding is always in equilibrium. 
However, Cherry and Adler demonstrate under a mild 
restriction that the fixed points for Eqs. ||5J) correspond 
to the possible fixed points of the full system. The mild 
restriction is that, under the full dynamics, if the concen- 
tration of one of the TFs is held fixed, the concentration 
of the other TF should settle on a unique value. This 
holds for our switch problem. 

Switching behaviour corresponds to the existence of 
two stable fixed points for Eqs. l|HJl, separated by a third 
fixed point which is unstable in one direction, like a 
transition state. Thus a test for switching behaviour 
is whether the dynamical system in Eqs. lO admits an 
unstable fixed point which is unstable in one direction. 
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From dynamical systems theory, this can be determined 
by considering the determinant of the stabihty matrix 



M 



fx ^ P-A fy 



9x 



9y - /^B 



(10) 



where — df/dx, etc. The required test is that 
det M < at the fixed point in question. At a fixed 
point, one can ehminate /Ia and ji^ between Eqs. ^ 
and ifTIUl . We find that detM = {'^/xy) x D where we 
have introduced the discriminant function 

D{x,y) = if ~ xf^){g - ygy) - {xg^)(ygy). (11) 

The sign of D mirrors the sign of det M so that ii D < 
at a fixed point, that fixed point is unstable in one 
direction. Moreover, for any point in the [x, y) plane, D 
has a definite value independent of the values of /2a and 
/2b- Thus in the plane of reduced concentration variables, 
the nature of a fixed point is determined solely by its 
position. In particular, the region D < Q contains all the 
unstable fixed points of the kind desired, and only those 
kind of fixed points. This region is bounded by the line 
£1 0. 

Normally one regards the parameters k^, etc, as given, 
and one has to solve Eqs. ijHJl for the position of any 
fixed points. However, if one knows a fixed point (x, y) 
in the plane of reduced concentration variables, the cor- 
responding reduced parameters are given by jlp^ = f /x 
and ji-Q — g/y, from Eqs. 0. The region 13 < in the 
(x, y) plane is therefore mapped to a region (a wedge) 
in the [jlp^,jlB) plane. Fig. |21 shows two examples. Note 
that the mapping from {x,y) to (/2a, /2b) is triple- valued 
within the wedge, since for parameters where switching 
behaviour occurs there are two stable fixed points in ad- 
dition to the unstable fixed point. 

An interesting analogy with a fluid-fluid demixing 
transition in a binary liquid mixture can be made at this 
point since both are examples of a cusp catastrophe 15]. 
In this analogy, jl^ and /2g^ correspond to the chemical 
potentials of the two components and the wedge D < to 
the region of spinodal instability. The cusp of the wedge 
in Fig.|2Ib) would correspond to the critical point. 

We now apply the above analysis to the switch models 
of the preceding section. Firstly, for n = m = 1 , one can 
prove that I? > for all four kinds of switch in Table 
This confirms the result of Cherry and Adler, namely 
that some form of co-operativity is required to make a 
working switch. 

Secondly, the discriminant for the totally co-operative 
case turns out to simplify to 13 = (1 -I- z)[l -f (m -t- n)z] 
where z = a;"y™. This is positive for all values of n and 
m. Thus switching behaviour cannot occur for the totally 
co-operative switch. 

For the remaining cases, we have analysed in detail the 
situation for dimerising switches, with n = m = 2. The 
details of this analysis are given in Appendix El Fig- 12] 
shows the regions in the both the {x,y) and {p,A, P-b) 
planes where switching behaviour occurs, for the general 




(a) 




FIG. 2: In mean field theory, switching behaviour is confined 
to a region in the {x, y) plane in (a), or equivalently to a wedge 
in the (/tA, Hb) plane in (b). Results are shown for dimerising 
general (solid line) and exclusive (dashed line) switches. Note 
that the bistable region is larger for the exclusive switch than 
for the general switch. 



and exclusive switches. Clearly there is a more extensive 
region of switching behaviour in the (/2a, /2b) plane for 
the exclusive case compared to the general case. Switch- 
ing behaviour is strongly suppressed for the partially co- 
operative cases; for example if OB2 is disallowed, switch- 
ing occurs only if /ja ^ 0.1 and /2b ^ 0.01 (lines not 
shown in Fig. Thus we conclude that, at least in 
mean field theory, restricting the set of operator states 
can have a marked effect on the possibility to form a 
genetic switch. 



IV. STOCHASTIC SIMULATIONS 

Whilst mean field theory indicates the regions in pa- 
rameter space where switching might occur, it has noth- 
ing to say about the switch stability and the effects of 
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TABLE II: Rate coefficients and equilibrium constants for a 
representative set of the reactions which define the baseline 
model. Here Vc is the cell volume and k is used as a unit of 
(inverse) time. 



reaction 


baseline 






O ^ + A 


k 








A ^ 




(0.2-0.8) A; 






2A ^ Aa 


fcf 


5kVc 






Aa ^ 2A 


kh 


5k 




= fcf/fcb = K 


O + Aa ^ OAa 


kon 


5kVc 






OAa ^ O + Aa 


fcoff 


k 







number fluctuations in systems with a finite size. To ad- 
dress these problems, we therefore turn to a stochastic 
simulation of the chemical reactions in Eqs. (Q. This 
is done via the 'Gillespie' algorithm jlS] , which is a ki- 
netic Monte Carlo scheme jl3| to generate trajectories in 
the space of concentration variables appropriate to the 
chemical master equation. 

Since these simulations are quite intensive, we focus 
again on the dimerising general and exclusive switches, 
and on the symmetry line = Mb = M ^^'^ = = 
k. We additionally have to supply values for all the rate 
coefficients in Eqs. To do this, we regard k « 0.1- 
1 s~"^ (the expression rate) as a unit of (inverse) time, and 
fi (the degradation rate) as the main control variable. As 
is the case with biological systems, real rate coefficients 
vary over quite a range of values. The baseline param- 
eter set is biologically motivated, with expression being 
a slow step and binding equilibrium biased in favour of 
bound states. Rate coefficients for our baseline model for 
a representative set of reactions are given in Table ITU For 
comparison, literature values for A-phage can be found in 
Arkin et al , in Aurell et al [zl la| , and in Bundschuh 
et al The rate coefficients for bimolecular reactions, 
such as the forward dimerisation and forward binding re- 
actions, have units of volume divided by time. Since our 
concentrations are expressed as the number of molecules 
in the cell volume, the natural units for these rate co- 
efficients are kVc where Vc ~ 2/im'^ ~ InM) is 
the cell volume. For the baseline parameters, the mean- 
field theory predicts bistability for fi/k < ^5/2 = 1.12 
for the general switch and fji/k < 2-^/5/3 = 1.49 for the 
exclusive switch. The implementation of the Gillespie al- 
gorithm with concentrations expressed as the number of 
molecules in the cell is straightforward. 

In the next section, we will first focus on the steady- 
state properties of both switches. In the subsequent sec- 
tion, we will address the dynamics of switching. 



A. 'Free-energy' landscapes 

We monitor the total number of molecules of A and B 
in the cell, TVa and TVb. These include the monomers. 
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FIG. 3: (a) Probabilities in {Na,Nb) plane constructed from 
trajectories of total duration kt m 8 x 10®, for general and 
exclusive switches, at /i = 0.4 k. Greyscale is logarithmic 
from P < 10~^ (white) to P > 0.04 (black), (b) Probabilities 
collapsed onto the Na — Nb line, plotted as a dimensionless 
'free energy' — log P; the maximum at A'^a ~ Nb corresponds 
to the mean-field transition state, (c) Probabilities along the 
diagonal line A^a = A'b; the minimum in this section now 
corresponds to the mean-field transition state. It is seen that 
the saddle-point of the exclusive switch is higher than that 
of the general switch. See also Table ITTTl for locations of the 
features in P{Na, Nb) and the barrier heights. 



the free dimers, and the bound dimcrs (see for example 
Eq. (|A15(l ). This is motivated primarily by the observa- 
tion by Bialek that the switch lifetime is likely to be an 
exponential function of the number of molecules involved 
in the switching process 'l9'|. Moreover, A'a and Ab are 
the relevant slow variables in the sense of Bundschuh et 
al 18]. 

To demonstrate the appearance of switching, we can 
construct the probability distribution P{Na,Nb) for 
states in the {Na,Nb) plane. This is done by generat- 
ing a sufficiently long Gillespie trajectory to obtain good 
coverage of all the interesting features of the probabil- 
ity distribution jT3|. Fig. 01 shows the main features for 
fj, = 0.4 A:, for both the dimerising general and exclusive 
switches with this direct method. 

We see that switching appears as a double maxi- 
mum in the probability distribution, and there is a 
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TABLE III: Locations of maxima and transition states for 
P{Na, Nb) for dimerising general (DGS) and exclusive (DES) 
switches ai fi = 0.4 (Fig-EJi compared to fixed points from 
mean field theory (MET, section ITTTl and appendix IXt. 









Maximum 


Transition state 








iVA 


Nb 


A^A = A^B 


DGS 




Gillespie 


15.0 ±0.5 





2 ± 1 


(7 = 


1) 


MET 


16.91 


0.09 


5.74 


DES 




Gillespie 


15.0 ±0.5 





4± 1 


(7 = 


0) 


MET 


16.04 


0.16 


3.15 



transition state at a low number of molecules for both 
TFs; we assume here that the transition state corre- 
sponds to the saddle point in the 'free-energy landscape' 
— \og[P{NA, Nb)]- Table UTTl contains the locations of 
the probability maxima and the transition state, com- 
paring the Gillespie results with the mean-field theory 
fixed point solutions of the preceeding section. The max- 
ima are in quite good agreement with the mean-field fixed 
points, although there is some difference in the location 
of the saddle-point. The first column in Table Hvl shows 
the probability of reaching states with Na — N-q deter- 
mined from FigO This probability is an order of magni- 
tude lower for the exclusive switch than for the general 
switch, which indicates that the exclusive switch is more 
stable than the general switch. 

Fig. 2] shows the mean number of molecules of the 
most-expressed transcription factor as a function of /i/fc 
in the switching regime, for the two kinds of switch. We 
define this via the time average 

]V = max(A^A,A^B). (12) 

Since the system spends most of the time near a sta- 
ble state, this average is dominated by the number of 
molecules in the stable states. As such, the time aver- 
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1 
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FIG. 4: Mean copy number of most-expressed TE, as a func- 
tion of the degradation rate, for both kinds of switch and 
baseline parameters. 'Gillespie' refers to the results of the 
stochastic simulations and 'MET' refers to those of the mean- 
field analysis described in section UTTl 




EIG. 5: Typical time series for TVa — Ab. Data is for the 
dimerising exclusive switch at jj, — OAk. 

age can be compared to the mean field stable fixed point. 
Further, in the stable steady states both switches should 
behave very similarly (see Eq.^. Fig.0]shows that while 
the average number of molecules A'^ is indeed nearly the 
same for the two switches, it also lies systematically a 
little above the mean field theory prediction. 

B. Rate constants 

We now turn to the switch dynamics. Fig.|Slshows typ- 
ical time series for A^a ^ for the dimerising exclusive 
switch at /i = 0.4 fc. The appearance of the two switch 
states where one of the TFs is strongly repressed com- 
pared to the other one is clearly seen. These correspond 
to the probability maxima of Fig.l^Ja). A switching event 
occurs when the roles of the two TFs flip spontaneously. 

In order to elucidate the stability of the switches, we 
make a connection with the reactive flux method that 
has been pioneered by Bennett 2Ql and Chandler 
and that is now widely used in the field of soft-condensed 
matter physics We first define an order parameter, 
q, that serves as our reaction coordinate and measures the 
progress of flipping the switch from one state to the other. 
In what follows, we will take q — Na — Ng. Furthermore, 
we define two characteristic functions that indicate in 
which state the system is in: 

hA{q) - 1, g < g* (13) 

= 0, q>q* (14) 

heiq) - 0, q<q* (15) 

= 1, q>q*. (16) 

It is natural to take for q* the value that corresponds 
to the top of the 'barrier' that separates the two sta- 
ble steady sates. We have thus chosen q* = for both 
switches (see Fig.|3)l. 

We now assume that we can make a separation of time 
scales [22. In particular, we assume that there exists 
a time t that is longer than the time itrans it takes for 
transient behavior to relax, but shorter than the charac- 
teristic time tj-xn for making a transition from one state 
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of the switch to the other. If this assumption holds, then 
we can coarse-grain the switch as a two-state system. If 
the flipping on the time scale trxn is a Poisson process, 
then the relaxation of the switch is given by the following 
expression 



M^M^^ = 7^(1 - eM-{kAB + kBA)t]) . (17) 
riA 

Here, the overbar denotes a time average; contributions 
to the average are weighted according to the stationary 
distribution of states. The expression on the left-hand 
side yields the probability that the switch is in state B 
at time t, given that initially it was in state A. If this 
time t is longer than itrans and if there are no correla- 
tions between successive switching events on this time 
scale, then this correlation function should be given by 
the macroscopic expression on the right-hand side of the 
above equation, where fc^s and ksA are the rate con- 
stants for flipping the switch in the forward and back- 
ward directions, respectively. The above relation thus 

only holds for t > ttrans- 

It is instructive to take the time derivative in Eq. 1171 

This yields 



^^^^^^^^^ = hBikAB + kBA)eM-ikAB + kBA)t]. (18) 

hA 

If we now consider times t smaller than ^ixn = I/I^ab + 
kBA), but larger than ttrans, and exploit the detailed bal- 
ance condition hB/hA = kAs/kBA, then we find 



hA{0)hB{t) 



= = kAB- 

Ha 



(19) 



This shows that the flux of trajectories from A to B 
should exhibit a plateau for times itrans < t < l/(kAB + 
kBA)- Indeed, the macroscopic switching rate is precisely 
given by the constant flux of trajectories in this regime. 

In general, however, it is useful to define a time depen- 
dent rate constant: 



kAB{t) 



hA{0)hB{t) 



hA 



m^m - - q*) 



oiq* - qit)) 



(20) 
(21) 



Siq-q*) 9(0)%(0) - g*)%(t) 



Oiq* -qit)) 



Siq~q*) 



= Poiq*) mOiqit) - q*) 
= Poiq*) Rit). 



(23) 
(24) 



Here, 9 is the Heaviside function. The overline with the 
asterix * denotes an average over states at the top of the 
barrier. 

It is seen that fc(<) = kAB it) is the product of two 
contributions. The first is Poiq*), which is given by, 



Poiq*) 



Piq* 



Eloo Piq) 



(25) 



Here Piq) is the steady-state probability of finding the 
system with a reaction coordinate of size q. Noting that 
\i q < q* the system is in the initial state, it is clear that 
Poiq*) is the probability of finding the system at the top 
of the barrier divided by the probability of finding it in 
the initial state. This quantity can be directly obtained 
from Fig. O and is shown in Table IIVI 

The second contribution to kit) is Rit), which gives 
the average flux of trajectories at the top of the barrier. 
As explained in [23|, the initial rate kit — > 0+) corre- 
sponds to the approximation for the rate constant in the 
transition state theory of chemical reactions: 



k 



■TST 



kit ^ 0+) ^ Poiq*) mOiq) 



(26) 



Transition-state theory assumes that all trajectories ini- 
tially heading from the top of the barrier towards the final 
state will indeed end up in the final state and all trajecto- 
ries initially heading towards the initial state will end up 
in the initial state. This assumption is only correct if no 
trajectories recross the barrier. In the present case, re- 
crossing turns out to be quite significant and, as a result, 
kit) will be significantly lower than fcxsT- It is conven- 
tional to express the reduction of fc(i) due to recrossings 
in terms of the transmission coefficient nit), defined as 



Kit) 



kit) _ Rit) 



i?(0^ 



(27) 



where -R(0+) — qiO)diq) . We note that, just as kABit) 
exhibits a plateau value for itrans <t < i^xn — ^/{kAB + 
kBA), i^it) and Rit) also reach a constant value on this 
time scale. We will simply refer to them as the transmis- 
sion coefficient k and kinetic prefactor R, respectively. 

For many rare event problems in soft-condensed mat- 
ter, it is not feasible to obtain the rate constant by cal- 
culating the correlation function (/iyi(0)/iB(t)) in a long, 
brute-force simulation. In order to alleviate the rare 
event problem, a widely used approach is to first calcu- 
late the free-energy barrier Poiq*) using umbrella sam- 
pling [23| and then compute the kinetic prefactor Rit) 
by shooting-off (molecular dynamics) trajectories from 
the top of the barrier ,2^ ,21j. However, the umbrella 
sampling technique in its conventional form relies on an 
energy-functional and is therefore only applicable to sys- 
tems that obey detailed balance. Biochemical networks 
are usually out of equilibrium and consequently lack de- 
tailed balance. In Appendix B, we show that the genetic 
toggle switches considered here also do not obey detailed 
balance. We can therefore not use the Bennett-Chandler 
approach [23, 0|. Indeed, we have resorted to long, 
brute-force simulations in order to calculate the flipping 
rates for the toggle switches. To be more precise we 
have obt ained the ra te constant by fitting the correlation 
function /iyi(0)/is(i)//iA to its macroscopic expression as 
given by Ea ll7l for the symmetric switches considered 
here, kAB — kBA — l/''"j where r is the life time of the 
stable steady state. In practice, to reduce noise in the 
correlation function, we excise a window of states around 
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FIG. 6: The cross correlation between the indicator func- 
tions hA = d{qA — q{t)) and hs — digit) — 9b) allows an 
accurate determination of the switch lifetime. Results are 
shown for the general and exclusive switches at /j, = 0.4 fc, 
with —QA — qs = 5. The lines are fits to Eq. in the form 
(/iA(0)/is(t)>/(/iA> = (1 - exp[-2f/r])/2. Table HVI contains 
the corresponding values of r. 



q* and define Ha = ^(<Za — qit)) and Hb = Oiqit) — qs) 
where qA < q* and qs > q* ■ Most of our calculations 
have been performed for —qA = qs = 5, but we find 
that, as expected, our results are insensitive to the pre- 
cise values chosen, provided that the system is well into 
the switching regime where there is a good separation of 
time scales. Fig. |B1 shows typical correlation functions. 

In Fig. Ul we show the switch lifetime r — l/kAs 
i— l/ksA for the symmetric switches) as a function of 
the mean value of the most-expressed TF for both kinds 
of switches, as /i varies. We see that r grows extremely 
rapidly with TV, which is the basic reason why extremely 
stable switches can be built with only a few hundred ex- 
pressed proteins. Our simulations cover 10 < < 30, 
but if we extrapolate our results to A'^ « 100, then 
kr w 10'' for the general switch but kr « 10^^ for the 
exclusive switch. In the latter case, this corresponds to 
lifetimes measured in tens of years. Such extremely long 
lifetimes have been reported for A-phage ||8| . Equally im- 
portant, Fig.[7|is a dramatic confirmation that the switch 
construction has a marked influence on stability. The ex- 
clusive switch lifetime grows much more rapidly with the 
mean copy number than the general switch lifetime. 

We find that r grows sub-exponentially with N and can 
be fit to a functional form suggested by analysis of a re- 
lated problem, namely that of switching between broken- 
symmetry phases in a driven diffusive model |24Ll25j |. The 
fit is to r = AlV"exp(6lV), where A = 50.4, a = 0.72, 
b = 0.097 for the general switch, and A = 7.32, a — 1.16, 
6 = 0.19 for the exclusive switch. This fit gives quanti- 
tative support to Bialek's conjecture that the switch life- 
time grows exponentially with the number of molecules 
involved in the switching process [19|, but additionally 
indicates that there is a logarithmic correction in N. 

Table HVl shows the barrier heights Poiq*) computed 
from the direct numerical simulations of PiNA, Nb). We 



TABLE IV: Switching kinetics. The first column is the prob- 
ability of reaching the top of the barrier from Fig. |3] the 
second column is the lifetime from Fig. |H| and the third col- 
umn is kinetic prefactor defined to be 7? = l/(Po (<?*)''"). The 
fourth column is the mean escape rate -R(0''') to A'^a > A'b 
from A^A = Nb as determined from additional analysis of the 
simulations, and the final column is the transmission coeffi- 
cient K = R/RiO+). A fi gure in brackets after a result is an 
estimate of the error in the final digit of that result. 



Po(g*)/10-^ 




R/k 


RiO+)/k 


K 


DGS 9(1) 


2.3(2) 


0.048(7) 


0.24(1) 


0.20(3) 


DES 0.92(1) 


8.0(5) 


0.14(1) 


0.98(4) 


0.14(1) 



can now obtain the kinetic prefactor R by dividing the 
rate constant kAB = 1/t by Poiq*)- R = kAB / Po{q*) 
(see also Eq. I24f) . Wc find that the underlying barrier- 
crossing rate R k, 0.05 — 0.15fc is a (small) fraction of 
k, which corresponds to the time scale for gene expres- 
sion (see Table In fact, for the general switch, the 
kinetic prefactor R is significantly lower than the slowest 
reaction in the system, which is the degradation reaction 
with ^ = 0.2 - 0.8fc (see also Table El). Since R can 
be interpreted as the rate at which the barrier is crossed 
(i.e., the flux of trajectories at the top of the barrier), it 
is clear that crossing the barrier typically involves a large 
number of protein synthesis and degradation steps. 

We have also directly computed from the Gillespie tra- 
jectories the distribution of escape times from the divid- 
ing surface A^a = -^b to states with A^a .^b- We 
flnd that for both kinds of switches, the escape time dis- 
tribution is quite well approximated by a Poisson dis- 
tribution with a characteristic lifetime r^. Since the 
switches are symmetric, this means we can define the 
escape rate to one side (for example to A^a > -^b) 
to be i?(0+) = 1/2t^. Estimates for i?(0+) from the 
Poisson-distribution fit are shown in the fourth column 



general 

exclusive jf' 




10 20 30 



N 



FIG. 7: Switch lifetime as a function of the mean number 
of the most-expressed TF. The exclusive switch becomes or- 
ders of magnitude more stable than the general switch at 
high numbers of the expressed TF. The lines are fits to 
T ~ N°'exp{bN) as discussed in the text. 
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in Table IIVI The ratio of the barrier crossing rate to 
the one-sided escape rate is the transmission coefficient, 
K = R/R{0+) 22]. As Table Irvl shows, the transmission 
coefficient for both kinds of switch is 10-20%. This shows 
that the barrier crossing is quite diffusive. In other words, 
a typical transition path between the two stable regions 
crosses and re-crosses the top of the barrier A'a = -^Vb sev- 
eral times before committing to one of the stable basins. 
This is in accord with our findings on the statistics of 
q = crossing times reported previously [inj |. It is also 
consistent with the observation that a typical transition 
path involves many protein synthesis and degradation re- 
actions. 



V. BEYOND THE BASELINE MODEL 

We now consider the effects of varying parameters 
away from the baseline set, and changing some other as- 
pects such as introducing a more detailed representation 
of transcription and translation. Switch stability is gen- 
erally unaffected by variation of the kinetic rates that do 
not alter the mean copy number of the TFs. What does 
have an effect, as we shall see, are phenomena which in- 
crease the 'shot noise' in expression. These include the 
generation of multiple copies of the TF from each mRNA 
transcript and stochastic fluctuations in the number of 
copies of the TF generated per mRNA transcript. 

In this section we give results for the dimerising ex- 
clusive switch. We have repeated all the simulations 
for the dimerising general switch and we find that the 
same trends are recovered. We also find that the ex- 
clusive switch is always markedly more stable than the 
corresponding general switch, in a manner represented 
by Fig. El 



A. Varying parameters away from baseline 

Fig. IHl shows that the switch lifetime is very insensi- 
tive to variations away from the baseline kinetic param- 
eters that do not affect the mean field steady states. If, 
however, we vary the kinetic parameters such that the 
mean field steady states are changed, then we do see a 
systematic change in the switch stability. For instance, 
the switch stability is enhanced if the binding equilib- 
rium is moved in favour of bound TF by increasing fcon 
(the rate of binding to DNA) . This is because the num- 
ber of molecules of each TF at the point where switch- 
ing just starts decreases as = kon/kos increases (see 
Eq. ljA16|l l and thus at fixed N we are effectively mov- 
ing deeper into the bistable region. This leads to a more 
stable switch. 




10 20 30 



N 

FIG. 8: Switch lifetime for the dimerising exclusive switch 
as some of the baseline kinetic parameters are varied. Data 
is shown for baseline model (circles; see Table ED, with the 
number of copies of the genome doubled and expression rate 
halved (squares), with kt and kh doubled (up triangles), with 
fcon and fcoff doubled (down triangles), with fcon = WkVc and, 
as before, fcoff = k (left triangles), and with fcon = 2kVc and, 
as before, fcofi — k (right triangles) ; see Table HTl for the mean- 
ing of the rate constants. The solid line is the fit to the base- 
line data from Fig. 

B. Effect of messenger RNA 

We have condensed the various steps in gene expression 
into a single 'expression reaction' in Eqs. Ijld|l . In reality, 
however, many steps are involved. We now analyse these 
steps in more detail. Recall that the genetic information 
is first transcribed into messenger RNA (mRNA), then 
translated into proteins by ribosomes. It is often the 
case that several copies of the protein can be generated 
from one mRNA transcript. In the simplest way, we can 
capture this by replacing Eqs. Ijld|l by a version in which 
r copies of each transcription factor are generated from 
each mRNA transcript, thus 

O I 0A„ rA X A, O I 0B,„ tb x B. (28) 

It is easy to see that nothing changes in mean field theory 
provided that rAkA and rBA:B are used in place of fcA and 
/cb. We have undertaken simulations for the case ta — 
tb = 2 and A:a = fcs = k/2, which in mean field theory 
is identical to our baseline. However, as Fig. |31 shows, 
the switch lifetime is decreased as a result of having TFs 
generated two at a time. The basic reason is that the 
'shot noise' in the expression reactions has increased. 

Even if on average one TF is produced from each 
mRNA transcript, noise can arise from the fluctuations 
about this mean, as we now discuss. To capture this, 
we explicitly include the generation and degradation of 
mRNA for the two TFs in the list of chemical reactions. 
The expression reactions in Eq. I|ld|l are replaced by the 
following reactions: 

O I 0A„ Ma ^ 0, (k'A), (ma) (29a) 

O I 0B„ Mb ^ 0, ik{,), iii'^) (29b) 

Ma-^A, Mb-^B. ik'l),ik'^) (29c) 
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FIG. 9; Switch stability is degraded if multiple copies of the 
TF are produced in each expression reaction, in this case 
rA = rs = 2 in Eq. II28|I . In a model which includes ex- 
plicit mRNA production (as described by Eqs. 1291 ). then 
even if on average only one TF is produced per transcript, 
the stability is similarly degraded by fluctuations about the 
mean. 'Baseline' refers to the baseline model, as described in 
Table El 

In these, Ma and Mb are the mRNA species, and the 
reactions represent the generation and degradation of 
mRNA, and the translation of mRNA into TFs. We have 
introduced new rate coefficients for the translation, tran- 
scription and mRNA degradation steps. 

In a steady state in mean field theory, the rate of tran- 
scription balances the rate of mRNA degradation and the 
mean number of mRNA transcripts is fc'//i'. Each tran- 
script produces TFs at a rate fc", hence the overall rate 
of production of TF is x k" . We conclude that 

the equivalent expression rate is 

k = (30) 

Going beyond mean field theory, the statistics of this 
basic model for mRNA translation can be solved ex- 
actly I23. Translation occurs at a rate k" and 
degradation at a rate hence the probability per unit 
time of either a translation or a degradation event tak- 
ing place is k" + /i'. Given that an event has taken 
place, the probability that it was a translation event is 
p = k" /{k" + fj,') and the probability that it was a degra- 
dation event is q = 1 — P — ^J-' /{k" + fi'). Therefore, 
the probability that r copies of the protein are generated 
before the mRNA transcript is degraded is Pr = qp^ ■ 
From this, all the statistics can be calculated. For in- 
stance, the mean number of TFs produced per tran- 
script is (r) = ^'^Q'rPr = p/q ^ k" / ^' , which agrees 
with the expectation from Eq. H30() . Also (r^) ^ p/q^, 
so that the standard deviation divided by the mean is 
\fqjp — (r)~^/^. Thus there is considerable noise due 
to fluctuations in the number of proteins generated per 
transcript, and this noise is largest when the mean num- 
ber of TFs produced per transcript is small. 

To see the effect of this noise we have implemented the 
additional reactions in Eqs. H29|l in our reaction schemes 



and again determined the switch lifetimes by direct sim- 
ulation. Fig. 1^ shows the results for k"f^ — k'^ = hk'p^ = 
Sfcg = [i!p^ = /ig = 5k. For these parameter values 
^a^a/a^b — ^b'^b/Mb — 1 s-iid, according to Eq. 1(20)), 
such a system should be identical to our baseline. More- 
over since k'J^/^'j^ = fegZ/ig — 1, there is on average one 
TF produced per mRNA transcript, so the 'shot noise' 
associated with the mean number of TFs produced per 
transcript is the same as in the baseline model. However, 
the simulations clearly demonstrate that the additional 
noise due to fluctuations about this mean reduces the 
switch lifetime considerably. 

C. Expression as a multistep process 

As a second refinement to the model for gene transcrip- 
tion and translation, we now consider gene expression as 
a multistep process. We thus replace the single reaction 
O ^ O + A by a sequence of reaction steps 

O^Oi^Oa^Os^ .0„^0 + A. (31) 

If there is only one intermediate stage, Oi say, this can 
be used to model the formation of an 'open complex' 
[28l | . If there are multiple intermediate stages, this could 
represent the individual steps of the RNA polymerase 
that walks along the DNA, or those of the ribosomes 
that walk along the mRNA. In some of the more detailed 
models of gene expression, all these intermediate steps 
are captured in detail 6] . 

We can make a connection with the baseline model 
by computing the waiting time for the lumped reaction 
O ^ O -I- A. Since the waiting times for the individual 
steps in Eq. H31|) are independent statistical quantities, 
the waiting time for the whole sequence is the sum of 
the waiting times for the individual steps. In terms of 
reaction rates, 1/^ = l/^ij where k is the rate of the 
lumped reaction, and the ki are the rates of the interme- 
diate steps. 

It follows that the waiting time distribution for the 
lumped reaction is not a Poisson distribution. Indeed the 
central limit theorem indicates that the lumped reaction 
will tend to have a Gaussian distribution of waiting times, 
effectively converging on a (5-function for a very large 
number of steps. Thus the approximation which replaces 
Eq. by O ^ O -|- A amounts to replacing the true 
non-Poisson distribution of waiting times by a Poisson 
distribution with the same mean. 

We have tested the consequences of this assumption 
for two cases. In the first test, we have inserted a single 
intermediate stage in Eq. H31() to represent the formation 
of the open complex. The rates of the two steps were 
chosen to be {5/A)k and 5k, so that formation of the open 
complex is the slow step. With this choice, the rate for 
the lumped reaction is the same as the baseline model. In 
the second test, we have inserted four intermediate stages 
in the reaction scheme, so that there are five intermediate 
reaction steps between O and 0-l-A, to model for example 
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FIG. 11: A 'split' exclusive switch can be built out of distinct 
operators provided each TF blocks the other TF from binding, 
at each operator site. 
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FIG. 10: Compared to the baseline model (circles; see Ta- 
ble [HJ, there is practically no effect on switch stability when 
one activation step is included (right triangles) or several are 
included (left triangles) as in Eq. (1311 (see main text for de- 
tails). There is also little change if the operator is 'split' as 
in Fig. inl and Eqs. 1321 (diamonds). The line is the fit to the 
baseline data from Fig. |7| 



the progressive stages of transcription. We have chosen 
the rates of these intermediate steps all equal to 5fc so 
again the effective rate for the lumped reaction is the 
same as the baseline model. However, the sum of the 
variances for the individual steps is now five times smaller 
than the variance in the waiting time for the lumped 
reaction. 

Fig. 111)1 shows data for the dimerising exclusive switch 
with these modifications. It shows that the switch sta- 
bility is practically unaffected by the inclusion of the ad- 
ditional reaction steps. These results suggest that the 
precise waiting time statistics for these multistcp reac- 
tions are less important for the switch stability than the 
statistics of the number fluctuations (as studied in con- 
siderable detail in the preceeding section). This is per- 
haps not so surprising, since the activation process to flip 
the switch must proceed by multiple coincidences of the 
'right' expression or degradation events (see also discus- 
sion on kinetic prefactor R in section IIV B|l . and so the 
detailed waiting time statistics of the individual reaction 
steps are unimportant. 



D. Can the operator be split? 

The basic premise of the exclusive switch is that the 
binding of one TF excludes the binding of the other TF. 
Whilst it is natural to think of this in terms of a diverging 
pair of genes on opposite strands of the DNA as shown 
in Fig. ^b), it is also possible to achieve the same effect 
by having each TF prevent the binding of the other TF 
at separated operators on the DNA. An exclusive switch 
with such a 'split' arrangement is shown in Fig. ^2 This 
arrangement is interesting because for each gene, one TF 
regulates expression by acting as a repressor, but the 
other TF acts only indirectly by preventing binding of 



the repressor TF. 

We can set up a system of chemical reactions to model 
the split operator case as follows. Suppose that Oa and 
Ob denote the operator sites for the genes which code for 
A and B, respectively. Then the set of reactions is 

Oa + A„^OaA„, Oa + B,„^OaB„, (32a) 
Ob + A„^ObA„, Ob + B„^ObB„, (32b) 
Oa I OaA„ a, Ob I ObB„ ^ B, (32c) 

plus the multimerisation and degradation reactions from 
Eqs. QJ. We assume that the rates for all binding and 
unbinding reactions are fcon and kos, and those for ex- 
pression of A and B are A:a, and k^, respectively. This 
set of reactions describes the split analogue of the exclu- 
sive switch. In mean field theory, this model is identical 
to the standard case. We have simulated this model and 
indeed find that there is little alteration to the switch 
stability, as shown in Fig. ^| 

We conclude that it is the exclusive nature of the TF 
binding which makes the exclusive switch markedly more 
stable than the general case. This can be achieved by 
overlapping the operator sites of genes that are arranged 
in a divergent manner and (thus) on opposite strands 
of the DNA as in Fig. Q^b), or by the pure exclusive 
binding arrangements of Fig. ^2 It appears, however, 
that nature has taken advantage of the elegance of the 
former scenario: in E. coli there are a significantly large 
number of diverging gene pairs with overlapping operator 
sites 



VI. DISCUSSION 

A genetic switch is inherently stochastic, because of 
the molecular character of its components. Our simula- 
tions demonstrate that the stability of a genetic switch 
can be strongly influenced by its construction. If the 
competing transcription factors mutually exclude each 
other at the operator regions, then the switch stability is 
markedly enhanced. Mutual exclusion of competing reg- 
ulatory molecules can be obtained by overlapping oper- 
ons, a network motif that we have recently identifled in 
a statistical analysis of the gene regulatory network of E. 
coli 113. 

The basic conclusion that mutual exclusion of compet- 
ing regulatory molecules can strongly enhance the sta- 
bility of biochemical networks is robust. Nevertheless, 
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the switch stabiHty is influenced by phenomena that in- 
crease the noise in gene expression. Such phenomena in- 
clude the generation of multiple copies of a transcription 
factor from the same mRNA transcript, as weh as intrin- 
sic noise arising from the fluctuations in the numbers of 
proteins produced from a transcript. 

When expressed as a function of the number of 
molecules involved in switching, the switch stability is 
characterised by a well deflned lifetime t which grows 
exponentially with the mean number, N, of the tran- 
scription factors that are involved in the switching. In 
the regime accessible to direct numerical simulations, the 
growth law is well characterised by r '--^ iV"exp(57V), 
where and a and b are parameters. 

Whereas we have investigated a number of details that 
might affect switch stability, we have left out some con- 
siderations which might additionally be important. Fore- 
most amongst these is the influence of cell division and 
the cell cycle. Other effects such as fluctuations in the 
availability of RNA polymerase and ribosomes might also 
have an influence on switch stability. Another aspect 
that should be investigated is the response of the switch 
to a perturbation, in other words how one might toggle 
a switch by introducing a pulse of some kind. We leave 
all these questions to future work. 

The rapid increase of the switch stability with the 
number of expressed transcription factors presents a fun- 
damental limitation to the use of direct simulation to 
compute the switch lifetime. This provokes the ques- 
tion: is there a smarter way to compute the switch sta- 
bility for long-lived switches? Such simulation methods 
would generically be useful since the characterisation of 
rare events remains a key challenge in computational sys- 
tems biology [23. In the field of soft-condensed matter 
physics, numerical techniques, such as the reactive flux 
method [23, 0] and transition path sampling [s^ , have 
been developed that make it possible to simulate rare 
events such as crystal nucleation, protein folding and 
chemical reactions. Biochemical networks, however, dif- 
fer fundamentally from these problems. In the former 
problems, the stationary distribution of states is usu- 
ally known beforehand. Indeed, these systems typically 
obey detailed balance. In contrast, our genetic circuits, 
like most biochemical networks, do not satisfy detailed 
balance. This means that the stationary distribution 
of states is not known a priori. As a result, numerical 
techniques developed to tackle rare events in the field of 
soft-condensed matter physics, cannot straightforwardly 
be transposed to biochemical networks. Recently, we 
have developed a new numerical technique, called For- 
ward Flux Sampling |3lj |. that does not rely upon prior 
knowledge of the stationary distribution of states. The 
scheme is easily applicable to a wide variety of biochem- 
ical switches, such as that of bacteriophage lambda, and 
makes it possible to calculate rates of switching and to 
identify patwhays for switching. The technique should 
therefore lead to a better understanding of these rare, 
but important events in biology. 
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APPENDIX A: MEAN FIELD THEORY FOR 
DIMERISING SWITCHES 

Here are the explicit results of the mean field (chemical 
kinetics) analysis for the dimerising switches. Firstly con- 
sider the general switch. For this case, a = /3 = 7 = 1, 
and Eqs. (0 and © simplify to / = 1/(1 + y^) and 
g = 1/(1 + x'^). The discriminant equation D = solves 
to 



and requires x > ~ 0.577 for solutions. The cusp 

of the wedge in the (/^a, Ab) plane lies at the point p.^ = 
/iB = 1/2 and corresponds to the point a; = y = 1 on the 
line D = 0. 

For the exclusive switch (a = /? = 1, 7 = 0) Eqs. (O 
and © give f{x,y) = (l-|-x^)/(H-x^ + j/^) and g{x,y) = 
f{y, x). The discriminant equation D = now solves to 



2 _ ±\/x^ + 14x6 + 25a;4 - 24a;2 -x'^ -bx^ + 2 
^ " 2{x^ - 1) 

(A2) 

This requires x > ^(3^17 - 11)/V2 ~ 0.827 for solu- 
tions, at which point y = -^3 ~ 1.73. The curve for 
X > 1 corresponds to the positive sign; one can show 
that the curve for x < 1 can be generated by reflecting 
the curve for a; > 1 through the line x — y (the curve 
passes through x — y ~ 1). Also, y — > 1 as a; — *■ 00. 
The cusp of the wedge in the (majMb) plane lies at the 
point //A = Mb = 2/3, again corresponding to the point 
X = y — I on the line D = 0. 

Finally consider the partially co-operative switch. This 
has a much reduced range of switching behaviour and 
the corresponding lines have not been shown in Fig. |21 
Consider the case where OB2 is disallowed. In this case 
(a = 7 = 1, /3 = 0) Eqs. © and © give / = (1 
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x^)/(l + + x^y^) and g = 1/(1 + + a;^y^). The 
discriminant equation D = Q solves to 



where 



y 



2 {^' + ^f 



x'^{x^ — 5) 



(A3) 



and requires x > y/5 « 2.24 to have a solution. 
The cusp of the wedge in the (/iAiMe) plane corre- 
sponds to the point x = ^(13 +Vl29)/V2 w 3.49 and 
y = 7(51 + 13V129)/10 w 1.41, and lies at /^a = 
v/(491 - 43v/129)/16 « 0.101 and /ig = V(159V129 - 
18057)/192 « 0.0190. 

We now describe the mean field fixed points for 
dimerising exclusive and general switches along the sym- 
metry line /iA = fJ-B = fJ- and = k^, = k. We assume 
there is one copy of the genome in the cell. The reduced 
degradation rate is 



A* = 



For the general switch, Eqs. ||HJ) are 

1 1 

L + 1 + 

The stable fixed points are 

1 ± - V 



{x,y) 



2fi 



(A5) 



(A6) 



(A7) 



where x > y ii the positive sign is taken for x. We see 
that jl < 1/2 is required as found earlier (see Fig. 
The unstable fixed point is 



X = y = 



u — 12/i^/u 
6/1 



where 



= lOS/i^ + 12/iV81 + 12/12. 



For the exclusive switch, Eqs. (O are 
1 -I- . 1 + 



1 -|- a; 2 4- J/ 2 ^ ' 1 + + y^ 
The stable fixed points are 

[l-2/i2 + ^l + 4^2±A]i/2 



My- 



(a;,?;) 



2/i 



(A8) 



(A9) 



(AlO) 



(All) 



where x > y ii the positive sign is taken for x, and 

= 2 - 12/7'* + 2(1 - 2/^2) Vl + 4/22. (^12) 

One can check that > requires jl < 2/3, as found 
earlier. The unstable fixed point is 



x^y 



1 + -y + (1 - 6fl^)/v 
6/1 



(A13) 



= 1 + Ahjl^ + 3/i V12 + 2T3yS2 + 24^ (A14) 

An expression for the total number of molecules of A in 
the cell which applies to both the general (7 = 1) and 
exclusive (7 = 0) switch is given by 



= riA + 2 nA2 + 2 tiqAs 
xVr 2x^Vr 



2{x^ +7a;2j/2) 



VKhKd K], 1 + x^ +y^ + -fx^y^ 

(A15) 

The total number of molecules of B is given by the same 
expression with x and y interchanged. This expression, 
together with the results above for the locations of the 
fixed points, is used to fill in the mean field rows in Ta- 
ble lllll For both kinds of switch, a; = y = 1 at the small- 
est value of jl for which bistability occurs. Therefore at 
this point. 



A^A = iVB = 



Vr. 



2Vc 



. 1 + 7 
'3-I-7' 



(A16) 



We finally consider briefly the limit /i 0, which is the 
limit of high expression and low degradation rates. In 
this limit, the stable fixed points for both exclusive and 
general switches converge with a; ~ 1/// and y ^ fi (or 
vice versa) . However the behavior of the unstable fixed 
point is X — y and x ^ ji^^/^ for the general switch, and 
X = y, and x ^ l/(2/i) for the exclusive switch. 



APPENDIX B: LACK OF DETAILED BALANCE 

We prove that the chemical master equation for the set 
of reactions in Eqs. 1^ cannot satisfy detailed balance. 
Recall that detailed balance implies that 



M^(l-> 2) Ps(l) = 1^(2-^1)^8(2) 



(Bl) 



where '1' and '2' indicate any two states, W{1 2) and 
14^(2 1) are the forward and backward transition rates 
between the states, and Ps is the steady-state probability 
distribution ■ If detailed balance holds, Ps is unique 
and all other probability distributions move towards Ps- 
For many problems, Ps is known (eg, it is the Gibbs- 
Boltzmann distribution). The problem is then to find 
a set of W which satisfies detailed balance, thus guar- 
anteeing that the dynamics has the correct equilibrium 
distribution. This is the case for Metropolis Monte-Carlo 
schemes for example. For the Gillespie algorithm or the 
chemical master equation however, the W are prescribed, 
and we do not know necessarily know Ps. 

Without prior knowledge of Ps, it would seem dif- 
ficult to determine from W alone whether a system 
obeys detailed balance or not. However Mukamel de- 
scribes a useful test |33|- Consider any cycle of states 
{1,2,3, ...,A:,1}. A necessary and sufficient condition 



14 



OBa ^ O + B2 OB2 ^ O + B2 

t i i T 

OB2+A ^ O + B2 + A OB2 + A -> O + B2+A (A4) 

fcoff X fcA X fcon X > Ox fcoff X /iA X fcon = 

FIG. 12: The fact that the products of the transition rates are different for the two cycles proves that the corresponding master 
equation cannot obey detailed balance. In both cases, the starting point is the state OB2. 



for the existence of detailed balance is that 

W{l^2)W{2^i)...W{k^l) 

^W{l-^k)W{k^k~l)...W{2^l) 

for all such cycles. To prove that a system does not 
obey detailed balance, it suffices to find a single counter- 
example. 

For Eqs. 1^, we can consider the cycles shown in 
Fig. El These cycles are mutual inverses, yet the cor- 
responding products of transition rates are clearly differ- 



ent. This proves that the master equation corresponding 
to Eqs. ^ does not obey detailed balance. 

It is clear that the same argument can be deployed 
whenever the appearance of a molecule in the system 
depends on the state of binding of some other molecules. 
This is very common, for instance it applies whenever one 
has gene regulation by a transcription factor. Thus the 
absence of detailed balance would appear to be a generic 
feature of any realistic biochemical reaction network. 
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